///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/fqmc/linear_method.cpp
///
/// \brief   implementation file for the linear method optimizer
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<vector>
#include<list>
#include<string>
#include<numeric>
#include<cassert>
#include<algorithm>
#include<cmath>

#include<boost/format.hpp>
#include<boost/shared_ptr.hpp>

#include<formic/fqmc/linear_method.h>
#include<formic/wfn/configuration.h>
#include<formic/fqmc/vmc_walker.h>
#include<formic/fqmc/energy.h>
#include<formic/wfn/wfn.h>
#include<formic/ham/ham.h>
#include<formic/input/base.h>
#include<formic/mpi/interface.h>
#include<formic/numeric/numeric.h>
#include<formic/numeric/matrix.h>
#include<formic/random/random.h>
#include<formic/output/output.h>
#include<formic/exception.h>
#include<formic/lapack/interface.h>
#include<formic/davidson.h>
#include<formic/genminres.h>
#include<formic/vector/vector.h>
#include<formic/vector/matrix.h>
#include<formic/timing/timing.h>
#include<formic/fqmc/linear_method_classes.h>
#include<formic/fqmc/line_search.h>
#include<formic/fqmc/energy_acc.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  constructs the linear method updater
///
/// \param[in]      userinp        the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> formic::fqmc::LinMethodUpdater<S>::LinMethodUpdater(const formic::InputBase & userinp)
  : _h_shifts(4, 0.0),
    _step_size(1.0)
{

  _h_shifts.at(3) = userinp.get<double>("lm_h_diag_shift");

}

template formic::fqmc::LinMethodUpdater<double>::LinMethodUpdater(const formic::InputBase & userinp);
template formic::fqmc::LinMethodUpdater<std::complex<double> >::LinMethodUpdater(const formic::InputBase & userinp);

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  performs a linear method update on the supplied wave function by a three-together
///         correlated sampling of updates' energies
///
/// \param[in]      userinp        the user's input options
/// \param[in]      input_sampler  object specifying the initial sample to use
/// \param[in,out]  vf             a shared pointer to the value function
///                                On exit, the best update has been added to the value function
///
/// \return Whether or not the update was rejected.  Updates are rejected if they do not
///         lower the energy or if there are problems in solving the linear method eigensystem.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> bool formic::fqmc::LinMethodUpdater<S>::update_3_together(const formic::InputBase & userinp,
                                                                             const boost::shared_ptr<formic::fqmc::Sampler<S> > & input_sampler,
                                                                             const boost::shared_ptr<formic::CoeffObj<S> > & vf) {

  // initialize the sampler for the LM solver
  boost::shared_ptr<formic::fqmc::Sampler<S> > sampler(input_sampler);

  // set the first set of shifts that will be tried
  _h_shifts.at(0) = _h_shifts.at(3) / 4.0;
  _h_shifts.at(1) = _h_shifts.at(3);
  _h_shifts.at(2) = _h_shifts.at(3) * 4.0;

  // get a vector to hold the energy/variance improvements obtained for different shifts
  std::vector<double> improvements(3, 0.0);

  // get a vector to hold the statistical errors of the improvements
  std::vector<double> stat_errs(3, 0.0);

  // get a vector to tell us which shifts have been solved
  std::vector<int> solved_shifts(3, false);

  // get a vector to hold the wave functions obtained for different shifts
  std::vector<boost::shared_ptr<formic::CoeffObj<S> > > shift_wfns(3);
  for (int i = 0; i < 3; i++)
    shift_wfns.at(i) = vf->clone();

  // save a copy of the wavefunction's initial variables
  formic::Vector<S,int> init_vars(vf->n_der(), formic::zero(S()));
  vf->to_vector(&init_vars.at(0));

  // compute initial energy, standard deviation, and statistical error
  double init_sdev = 0.0;
  double init_serr = 0.0;
  const double init_energy = formic::fqmc::energy(userinp, sampler, vf, init_sdev, init_serr, true);

  // create the eigensolver
  formic::fqmc::DavidsonLM<S> dsolver(userinp, vf, init_energy, init_sdev);

  // get a vector to hold the eigenvector
  formic::Vector<S,int> evec(1+vf->n_der(), formic::zero(S()));

  // set flag telling us if the solver's krylov subspace needs to be reset
  bool solver_needs_to_reset = true;

  // set a flag to tell if we have a set of three good shifts
  bool shift_success = false;

  // keep trying central shifts until one works
  for (int si = 0; si < 20 && !shift_success; si++) {

    // if the previous attempt failed, draw a new sample and use it for the davidson solver
    if ( si > 0 ) {

      // create a sampler that will take a new random sample
      boost::shared_ptr<formic::fqmc::Sampler<S> > temp_sampler( new formic::fqmc::Sampler<S>(userinp,
                                                                                              userinp.get<int>("energy_burn_len"),
                                                                                              userinp.get<int>("energy_samp_len"),
                                                                                              formic::Configuration(userinp),
                                                                                              formic::LaggedFibonacci<>(),
                                                                                              vf ) );
      sampler = temp_sampler;

      // reset the davidson solver so that it uses the new sample
      dsolver.reset(userinp, sampler, vf, _h_shifts.at(0), true); // last argument forces matrices to rebuild for new sample
      solver_needs_to_reset = false;

      // make sure the solution gets updated for each shift
      std::fill(solved_shifts.begin(), solved_shifts.end(), false);

    }

    // get the solution for each shift
    for (int p = 0; p < 3; p++) {

      // if we've already solved this shift, skip it
      if ( solved_shifts.at(p) )
        continue;

      // if necessary, reset the eigensolver and the eigenvector guess
      if ( solver_needs_to_reset || userinp.get<bool>("lm_restart_krylov") || dsolver.min_S_eval() < userinp.get<double>("lm_min_S_eval") ) {

        // reset eigensolver (for davidson this also empties the krylov subspace)
        dsolver.reset(userinp, sampler, vf, _h_shifts.at(p));
        solver_needs_to_reset = false;

        // for davidson, set initial guess to the current wave function
        evec = formic::zero(S());
        evec.at(0) = formic::unity(S());
        dsolver.add_krylov_vector(userinp, evec);

      }

      // print which shift's eigensystem is being solved
      if (formic::mpi::rank() == 0)
        formic::of << boost::format("solving linear method eigensystem for shift %.4e") % _h_shifts.at(p) << std::endl << std::endl;

      // update the eigensolver for the current shift
      dsolver.update_ham_shift(userinp, _h_shifts.at(p));

      // compute the lowest eigenvector
      double davidson_eval = 0.0;
      bool good_solve = dsolver.solve(userinp, davidson_eval, evec);

      // set the bad solve flag if the eigensolver found an eigenvector that is not dominated by the current state
      const double max_update_abs_value = std::abs(std::accumulate(evec.begin()+1, evec.end(), formic::zero(S()), formic::LargestAbsAccumulator<S>()));
      const double max_allowable_update_abs_value = userinp.get<double>("lm_max_update_abs");
      good_solve = ( good_solve && max_update_abs_value < max_allowable_update_abs_value );

      // print the largest weight on a derivative vector
      if ( formic::mpi::rank() == 0 )
        formic::of << boost::format("largest magnitude derivative vector weight for shift %.4e is %.6e") % _h_shifts.at(p) % max_update_abs_value << std::endl << std::endl;

      // send the flag telling if the solve was good to all processes
      formic::mpi::bcast(good_solve);

      // record whether we've solved the shift successfully
      solved_shifts.at(p) = good_solve;

      // if the solve succeeded, process the update
      if ( good_solve ) {

        // "normalize" the update and convert it into wave function variables
        if ( formic::mpi::rank() == 0 ) {
          evec /= 1.0 * evec[0];
          formic::xaxpy(vf->n_der(), formic::unity(S()), &init_vars.at(0), 1, &evec.at(1), 1);
        }

        // save the wave function variables resulting from the update
        formic::mpi::bcast(evec);
        shift_wfns.at(p)->from_vector(&evec.at(1));

      // otherwise, we need to reset the solver to clear the krylov subspace
      } else {
        solver_needs_to_reset = true;
      }

    }

    // if all solves were successful, compute energies with correlated sampling
    if ( solved_shifts.at(0) && solved_shifts.at(1) && solved_shifts.at(2) ) {

      // get energy accumulators for updated wave functions
      std::vector<boost::shared_ptr<formic::Accumulator<S> > > en_accs(3);
      //{
      //  boost::shared_ptr<formic::Accumulator<S> > eacc( new formic::fqmc::EnAcc<S>(shift_wfns.at(1),
      //                                                                              vf,
      //                                                                              userinp.get<int>("energy_ac_length"),
      //                                                                              userinp.get<bool>("exact_sampling"),
      //                                                                              userinp.get<bool>("print_config_info")) );
      //  en_accs.at(3) = eacc;
      //}
      for (int p = 0; p < 3; p++) {
        boost::shared_ptr<formic::Accumulator<S> > eacc( new formic::fqmc::EnAcc<S>(shift_wfns.at(1),
                                                                                    shift_wfns.at(p),
                                                                                    userinp.get<int>("energy_ac_length"),
                                                                                    userinp.get<bool>("exact_sampling"),
                                                                                    userinp.get<bool>("print_config_info")) );
        en_accs.at(p) = eacc;
      }

      // get energy accumulator for initial wave function
      boost::shared_ptr<formic::Accumulator<S> > eacc_init( new formic::fqmc::EnAcc<S>(shift_wfns.at(1),
                                                                                       vf,
                                                                                       userinp.get<int>("energy_ac_length"),
                                                                                       userinp.get<bool>("exact_sampling"),
                                                                                       userinp.get<bool>("print_config_info")) );

      // create sampler
      formic::fqmc::Sampler<S> shift_sampler(userinp,
                                             userinp.get<int>("energy_burn_len"),
                                             userinp.get<int>("energy_samp_len"),
                                             formic::Configuration(userinp),
                                             formic::LaggedFibonacci<>(),
                                             shift_wfns.at(1) );

      // sample the energies
      shift_sampler.simple_sample() << shift_wfns.at(0) << shift_wfns.at(2) << en_accs.at(0) << en_accs.at(1) << en_accs.at(2) << vf << eacc_init;
      //shift_sampler.simple_sample() << vf << shift_wfns.at(0) << shift_wfns.at(2) << en_accs.at(0) << en_accs.at(1) << en_accs.at(2) << en_accs.at(3);

      // get the statistical errors for each energy
      for (int p = 0; p < 3; p++)
        stat_errs.at(p) = ((formic::fqmc::EnAcc<S>*)(en_accs.at(p).get()))->serr();

      // check if the sampled energies' statistical errors are reasonable
      for (int p = 0; p < 3; p++)
       solved_shifts.at(p) = ( solved_shifts.at(p) && stat_errs.at(p) < userinp.get<double>("lm_stat_err_tol") );

      // get energies and check that at least one of the shifts' energies is reasonably low
      bool all_energies_too_high = true;
      for (int p = 0; p < 3; p++) {
        improvements.at(p) = ((formic::fqmc::EnAcc<S>*)(en_accs.at(p).get()))->energy();
        //if ( improvements.at(p) < init_energy + init_serr )
        if ( improvements.at(p) < ((formic::fqmc::EnAcc<S>*)(eacc_init.get()))->energy() + ((formic::fqmc::EnAcc<S>*)(eacc_init.get()))->serr() )
          all_energies_too_high = false;
      }
      formic::mpi::bcast(all_energies_too_high);

      // print results for the different shifts
      if ( formic::mpi::rank() == 0 ) {
        formic::of << std::endl;
        formic::of << boost::format("%20s      %20s      %20s") % "shift       " % "energy       " % "stat error     " << std::endl;
        formic::of << boost::format("--------------------      --------------------      --------------------") << std::endl;
        for (int p = 0; p < 3; p++)
          formic::of << boost::format("%20.6e      %20.12f      %20.6f") % _h_shifts.at(p) % improvements.at(p) % stat_errs.at(p) << std::endl;
        formic::of << boost::format("%20s      %20.12f      %20.6f")
                        % "initial values"
                        % ((formic::fqmc::EnAcc<S>*)(eacc_init.get()))->energy() // % init_energy
                        % ((formic::fqmc::EnAcc<S>*)(eacc_init.get()))->serr() // % init_serr
                   << std::endl;
        formic::of << std::endl;
      }

      // if no energies are reasonably low, move to higher shifts
      if ( all_energies_too_high ) {
        std::fill(solved_shifts.begin(), solved_shifts.end(), false);
        if ( formic::mpi::rank() == 0 )
          formic::of << boost::format("All shifts' energies were too high.  Will attempt larger shifts.") << std::endl << std::endl;
      }

    }

    // if any of the solvers failed, drop the lowest shift and try again with larger shifts
    if ( !solved_shifts.at(0) || !solved_shifts.at(1) || !solved_shifts.at(2) ) {
      for (int p = 0; p < 2; p++) {
        std::swap(_h_shifts.at(p), _h_shifts.at(p+1));
        std::swap(solved_shifts.at(p), solved_shifts.at(p+1));
        std::swap(shift_wfns.at(p), shift_wfns.at(p+1));
      }
      _h_shifts.at(2) = _h_shifts.at(1) * 2.0;
      solved_shifts.at(2) = false;
      shift_success = false;

    // otherwise, we are done moving the shifts
    } else {
      shift_success = true;
    }

  }

  // if we have found three successful shifts, interpolate to find the best shift
  int best_i = -1;
  if ( shift_success && formic::mpi::rank() == 0 ) {
    const double energy_thresh = 1.0e-9;
    if        ( improvements.at(0) < improvements.at(1) && improvements.at(0) < improvements.at(2) ) {
      best_i = 0;
    } else if ( improvements.at(2) < improvements.at(0) && improvements.at(2) < improvements.at(1) ) {
      best_i = 2;
    } else if ( std::abs( improvements.at(0) - improvements.at(1) ) < energy_thresh && std::abs( improvements.at(1) - improvements.at(2) ) < energy_thresh ) {
      best_i = 1;
    } else {
      std::vector<double> log_shifts(3, 0.0);
      for (int i = 0; i < log_shifts.size(); i++)
        log_shifts.at(i) = std::log(_h_shifts.at(i));
      _h_shifts.at(3) = std::exp( formic::porabola_min_max(&log_shifts.at(0), &improvements.at(0)) );
      if ( _h_shifts.at(3) > 1.1 * _h_shifts.at(2) )
        throw formic::Exception("unexpectedly large shift of %.4e found in the linear regression during the linear method") % _h_shifts.at(0);
      if ( _h_shifts.at(3) < 0.9 * _h_shifts.at(0) )
        throw formic::Exception("unexpectedly small shift of %.4e found in the linear regression during the linear method") % _h_shifts.at(0);
    }
    if ( best_i >= 0)
      _h_shifts.at(3) = _h_shifts.at(best_i);
    _h_shifts.at(3) = std::max(_h_shifts.at(3), userinp.get<double>("lm_min_h_diag_shift"));
  }

  // send shift and shift index to all processes
  if ( shift_success ) {
    formic::mpi::bcast(_h_shifts.at(3));
    formic::mpi::bcast(best_i);
  }

  // if we used parabolic interpolation, get the update for the resulting shift
  if ( shift_success && best_i < 0 ) {

    // if necessary, reset the eigensolver and the eigenvector guess
    if ( solver_needs_to_reset || userinp.get<bool>("lm_restart_krylov") || dsolver.min_S_eval() < userinp.get<double>("lm_min_S_eval") ) {

      // reset eigensolver (for davidson this also empties the krylov subspace)
      dsolver.reset(userinp, sampler, vf, _h_shifts.at(3));
      solver_needs_to_reset = false;

      // for davidson, set initial guess to the current wave function
      evec = formic::zero(S());
      evec.at(0) = formic::unity(S());
      dsolver.add_krylov_vector(userinp, evec);

    }

    // update the eigensolver for the current shift
    dsolver.update_ham_shift(userinp, _h_shifts.at(3));

    // compute the lowest eigenvector
    double davidson_eval = 0.0;
    bool good_solve = dsolver.solve(userinp, davidson_eval, evec);

    // "normalize" the update and convert it into wave function variables
    if ( formic::mpi::rank() == 0 && good_solve ) {
      evec /= 1.0 * evec[0];
      formic::xaxpy(vf->n_der(), formic::unity(S()), &init_vars.at(0), 1, &evec.at(1), 1);
    }

    // if the solve failed, we will reject the update
    if ( !good_solve ) {
      shift_success = false;
      if (formic::mpi::rank() == 0)
        formic::of << "linear method failed to solve the eigensystem for the interpolated shift and so will reject the update" << std::endl << std::endl;
    }

  // otherwise, get the wave function variables from the previously solved shift
  } else if ( shift_success && formic::mpi::rank() == 0 ) {
    shift_wfns.at(best_i)->to_vector(&evec.at(1));
  }

  // if we succeeded, update the value function
  bool rejected;
  if ( shift_success ) {

    rejected = false;
    formic::mpi::bcast(evec);
    vf->from_vector(&evec.at(1));

  // otherwise, reject the update
  } else {

    rejected = true;
    if (formic::mpi::rank() == 0)
      formic::of << "linear method is rejecting the update" << std::endl << std::endl;
    _h_shifts.at(3) = *std::max_element(_h_shifts.begin(), _h_shifts.end());
    formic::mpi::bcast(_h_shifts.at(3));

  }

  // print the suggested next shift
  if (formic::mpi::rank() == 0)
    formic::of << boost::format("suggested next shift is %.6e") % _h_shifts.at(3) << std::endl << std::endl;

  // return whether or not the update was rejected
  return rejected;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  performs a linear method update on the supplied wave function
///
/// \param[in]      userinp        the user's input options
/// \param[in]      sampler        object that controls sampling
/// \param[in,out]  vf             a shared pointer to the value function
///                                On exit, the best update has been added to the value function
///
/// \return Whether or not the update was rejected.  Updates are rejected if they do not
///         lower the energy.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> bool formic::fqmc::LinMethodUpdater<S>::update(const formic::InputBase & userinp,
                                                                  boost::shared_ptr<formic::fqmc::Sampler<S> > & sampler,
                                                                  const boost::shared_ptr<formic::CoeffObj<S> > & vf) {

  // if requested, find update by correlated sampling of three shifts' energies together
  if ( userinp.get<bool>("lm_three_together") )
    return this->update_3_together(userinp, sampler, vf);

  // set the shifts that will be tried
  _h_shifts.at(0) = _h_shifts.at(3) * 10.0;
  _h_shifts.at(1) = _h_shifts.at(3);
  _h_shifts.at(2) = _h_shifts.at(3) / 10.0;

  // get a vector to hold the energy/variance improvements obtained for different shifts
  std::vector<double> shift_improvements;

  //// get a vector to hold the wave functions obtained for different shifts
  //std::vector<std::vector<S> > shift_wfns(_h_shifts.size());
  //for (int i = 0; i < shift_wfns.size(); i++)
  //  shift_wfns.at(i).assign(vf->n_der(), formic::zero(S()));

  // save a copy of the wavefunction's initial variables
  formic::Vector<S,int> init_vars(vf->n_der(), formic::zero(S()));
  vf->to_vector(&init_vars.at(0));

  // compute energy and standard deviation of the initial wavefunction for this iteration
  double init_sdev = 0.0;
  double init_serr = 0.0;
  const double init_energy = formic::fqmc::energy(userinp, sampler, vf, init_sdev, init_serr);

  // create the eigensolver
  formic::fqmc::DavidsonLM<S> dsolver(userinp, vf, init_energy, init_sdev);

  // get a vector to hold the eigenvector
  formic::Vector<S,int> evec(1+vf->n_der(), formic::zero(S()));

//  // set the initial guess equal to the initial wavefunction
//  std::vector<S> evec(1+vf->n_der(), formic::zero(S()));
//  evec.at(0) = formic::unity(S());
//
//  if ( userinp.get<bool>("lm_do_davidson") && !userinp.get<bool>("lm_restart_krylov") ) {
//
//    // reset the davidson solver (this also empties the krylov subspace)
//    dsolver.reset(userinp, gf, vf, start_config, start_rng, _h_shifts.at(0));
//
//    // add the initial guess to the krylov subspace
//    dsolver.add_krylov_vector(userinp, evec);
//
//    //// add previous updates to the krylov subspace
//    //for (typename std::list<std::vector<S> >::iterator it = prev_updates.begin(); it != prev_updates.end(); it++)
//    //  dsolver.add_krylov_vector(userinp, *it);
//
//  } else if ( !userinp.get<bool>("lm_do_davidson") ) {
//
//    // reset the davidson solver
//    dsolver.reset(userinp, gf, vf, start_config, start_rng, _h_shifts.at(0));
//
//  }

  // get flag to tell if any of the shifts gave a bad solve
  bool any_bad_solves = false;

  // get flag to tell if the most recent shift gave a bad solve
  bool previous_solve_was_bad = false;

  // get variables to hold the best energy/variance improvement, best shift, and best variables
  double best_improvement = 1.0;
  double best_shift = _h_shifts.at(1);
  formic::Vector<S,int> best_vars = init_vars.clone();

  // compute a wavefunction update and its energy for each of the shifts
  for (int si = 0; si < _h_shifts.size() && si < 20; si++) {

    // make room for this shift's improvement
    shift_improvements.push_back(1.0);

    // print which shift is being tried
    if (formic::mpi::rank() == 0)
      formic::of << boost::format("trying a shift of %.4e") % _h_shifts.at(si) << std::endl << std::endl;

    // restore the wavefunction
    vf->from_vector(&init_vars.at(0));

    // if necessary, reset the eigensolver and the eigenvector guess
    if ( si == 0 || userinp.get<bool>("lm_restart_krylov") || dsolver.min_S_eval() < userinp.get<double>("lm_min_S_eval") || previous_solve_was_bad ) {

      // reset eigensolver (for davidson this also empties the krylov subspace)
      dsolver.reset(userinp, sampler, vf, _h_shifts.at(si));

      // for davidson, set initial guess to the current wave function
      evec = formic::zero(S());
      evec.at(0) = formic::unity(S());
      dsolver.add_krylov_vector(userinp, evec);

    }

    // update the eigensolver for the new shift
    dsolver.update_ham_shift(userinp, _h_shifts.at(si));

    // compute the lowest eigenvector
    double davidson_eval = 0.0;
    bool good_solve = dsolver.solve(userinp, davidson_eval, evec);

    // "normalize" the update
    if ( formic::mpi::rank() == 0 )
      evec /= evec[0];

    // send the update to all processes
    formic::mpi::bcast(evec);

    //// print the update
    //if ( formic::mpi::rank() == 0 ) {
    //  formic::Matrix<S,int> udm(evec.size(), 1);
    //  udm <<= &evec.at(0);
    //  std::sort(udm.begin()+1, udm.end(), formic::AscendingNormSorter<S>());
    //  formic::of << udm.print("%20.12f", "sorted update vector") << std::endl;
    //}

    // print the largest weight on a derivative vector
    if ( formic::mpi::rank() == 0 )
      formic::of << boost::format("largest magnitude derivative vector weight = %.6e")
                    % std::abs(std::accumulate(evec.begin()+1, evec.end(), formic::zero(S()), formic::LargestAbsAccumulator<S>()))
                 << std::endl
                 << std::endl;

    // get the variance and energy weightings for the optimization
    const double vw = userinp.get<double>("lm_var_weight");
    const double ew = 1.0 - vw;

    // set the bad solve flag if the eigensolver found an eigenvector that is not dominated by the current state
    const double max_update_abs_value = std::abs(std::accumulate(evec.begin()+1, evec.end(), formic::zero(S()), formic::LargestAbsAccumulator<S>()));
    const double max_allowable_update_abs_value = userinp.get<double>("lm_max_update_abs");
    good_solve = ( good_solve && max_update_abs_value < max_allowable_update_abs_value );

    // set the bad solve flag if the eigensolver found an eigenvalue that was too far from the current energy
    good_solve = ( good_solve && std::abs( davidson_eval - ew * init_energy - vw * init_sdev * init_sdev ) < userinp.get<double>("lm_max_e_change") );

    // set the bad solve flag if the statistical error of the update is too large
    formic::Vector<S,int> err_check_save(vf->n_der(), formic::zero(S()));
    vf->to_vector(&err_check_save.at(0));
    vf->update_variables(&evec.at(1),  formic::unity(S()) * _step_size);
    double shift_sdev = 1.0e100;
    double shift_serr = 1.0e100;
    double shift_energy = formic::fqmc::energy(userinp, sampler, vf, shift_sdev, shift_serr);
    vf->from_vector(&err_check_save.at(0));
    good_solve = ( good_solve && shift_serr < userinp.get<double>("lm_stat_err_tol") );

    //// if we are not doing a line search, then the solve is bad if the resulting standard deviation is too large
    //if ( !userinp.get<bool>("lm_do_line_search") )
    //  good_solve = ( good_solve && shift_sdev > userinp.get<double>("update_max_sdev") );

    // send the solver status flag to all processes
    formic::mpi::bcast(good_solve);

    // if the solver succeeded, use the eigenvector to compute an update
    if (good_solve) {

      //// print the eigenvector
      //if ( formic::mpi::rank() == 0 ) {
      //  for (int i = 0; i < evec.size(); i++)
      //    formic::of << boost::format("%20.12f") % evec.at(i) << std::endl;
      //  formic::of << std::endl;
      //}

      // if requested, perform a line search along the update direction
      if ( userinp.get<bool>("lm_do_line_search") ) {

        // do line search
        double step_arg = _step_size;
        formic::fqmc::line_search(userinp, &evec.at(1), userinp.get<double>("update_max_sdev"), sampler, vf, step_arg);

        // re-compute new energy and standard deviation using correlated sampling
        if (formic::mpi::rank() == 0)
          formic::of << boost::format("computing energy of updated wavefunction for the %.4e shift...") % _h_shifts.at(si) << std::endl << std::endl;
        shift_sdev = 1.0e100;
        shift_serr = 1.0e100;
        shift_energy = formic::fqmc::energy(userinp, sampler, vf, shift_sdev, shift_serr);

      // otherwise, assume that a unit step length is good
      } else {
        vf->update_variables(&evec.at(1),  formic::unity(S()) * _step_size);
      }

      // compute the improvement over the initial wave function
      shift_improvements.at(si) = ew * ( shift_energy - init_energy ) + vw * ( shift_sdev * shift_sdev - init_sdev * init_sdev );

      // if the standard deviation was too large, we act as if the solver was bad (this will force larger shifts to be tried)
      good_solve = ( good_solve && shift_sdev < userinp.get<double>("update_max_sdev") );

      // send the solver status flag to all processes
      formic::mpi::bcast(good_solve);

      // if this was the best shift so far, save it, its improvement, and its variables
      if ( good_solve && shift_improvements.at(si) < best_improvement ) {
        best_shift = _h_shifts.at(si);
        best_improvement = shift_improvements.at(si);
        vf->to_vector(&best_vars.at(0));
      }

    }
    
    // remember that this shift produced a healthy solution
    if ( good_solve )
      previous_solve_was_bad = false;

    // if the solver failed or the sdev was too large, abandon the current shift and try larger shifts until the problem goes away
    if (!good_solve) {

      // remember that there was a bad solve
      any_bad_solves = true;

      // remember that the most recent solve was bad
      previous_solve_was_bad = true;

      // set the next shift to be larger
      const double larger_shift = *std::max_element(_h_shifts.begin(), _h_shifts.begin()+si+1) * 10.0;
      if ( si+1 == _h_shifts.size() )
        _h_shifts.push_back(larger_shift);
      else if ( si+1 < _h_shifts.size() )
        _h_shifts.at(si+1) = larger_shift;
      else
        throw formic::Exception("_h_shifts vector had an unexpectedly small length of %i") % _h_shifts.size();
      formic::mpi::bcast(_h_shifts);

      // set the current improvement to be very bad
      shift_improvements.at(si) = 1.0e100;

      // declare what we are doing
      if (formic::mpi::rank() == 0)
        formic::of << boost::format("!!! Abandoning current bad shift.  Larger shifts will be tried until one works.") << std::endl << std::endl;

    }

    // send the improvement to all processes
    formic::mpi::bcast(shift_improvements.at(si));

    // if we have done three shifts and they were all good, make a guess for the best value of the fourth shift
    if ( si == 2 && !any_bad_solves && formic::mpi::rank() == 0 ) {
      const double energy_thresh = 1.0e-9;
      if        ( shift_improvements.at(0) < shift_improvements.at(1) && shift_improvements.at(0) < shift_improvements.at(2) ) {
        _h_shifts.at(3) = _h_shifts.at(0) * 3.0;
      } else if ( shift_improvements.at(2) < shift_improvements.at(0) && shift_improvements.at(2) < shift_improvements.at(1) ) {
        _h_shifts.at(3) = _h_shifts.at(2) / 3.0;
      } else if ( std::abs( shift_improvements.at(0) - shift_improvements.at(1) ) < energy_thresh && std::abs( shift_improvements.at(1) - shift_improvements.at(2) ) < energy_thresh ) {
        _h_shifts.at(3) = _h_shifts.at(1);
      } else {
        std::vector<double> log_shifts(3, 0.0);
        for (int i = 0; i < log_shifts.size(); i++)
          log_shifts.at(i) = std::log(_h_shifts.at(i));
        _h_shifts.at(3) = std::exp( formic::porabola_min_max(&log_shifts.at(0), &shift_improvements.at(0)) );
        if ( _h_shifts.at(3) > 1.1 * _h_shifts.at(0) )
          throw formic::Exception("unexpectedly large shift of %.4e found in the linear regression during the linear method") % _h_shifts.at(0);
        if ( _h_shifts.at(3) < 0.9 * _h_shifts.at(2) )
          throw formic::Exception("unexpectedly small shift of %.4e found in the linear regression during the linear method") % _h_shifts.at(0);
      }
      _h_shifts.at(3) = std::max(_h_shifts.at(3), userinp.get<double>("lm_min_h_diag_shift"));
    }

    // broadcast the shift list to all processes
    formic::mpi::bcast(_h_shifts);

  }

//  // save the update for use in seeding future krylov subspaces
//  if (formic::mpi::rank() == 0)
//    prev_updates.push_back(evec);
//  else
//    prev_updates.push_back(std::vector<S>(1));
//  prev_updates.rbegin()->at(0) = formic::zero(S());
//  if ( prev_updates.size() > max_prev_updates )
//    prev_updates.pop_front();

  bool rejected = false;
  double next_shift = 1.0;

  // if we found a shift that lowered the target function, update the wave function accordingly and check its energy on a fresh sample
  bool fresh_energy_acceptable = false;
  if ( best_improvement < 0.0 ) {
    fresh_energy_acceptable = true;
    vf->from_vector(&best_vars.at(0));
    //if (formic::mpi::rank() == 0)
    //  formic::of << "checking update's energy on a fresh sample" << std::endl << std::endl;
    //const double fresh_energy = formic::fqmc::energy(userinp, vf);
    //fresh_energy_acceptable = ( fresh_energy < init_energy + 2.0 * init_serr );
  }

  // if the fresh energy is ok, save the shift for next time
  if ( fresh_energy_acceptable ) {

    next_shift = best_shift;

  // otherwise, reject the update
  } else {
    rejected = true;
    if (formic::mpi::rank() == 0)
      formic::of << "linear method failed to lower the energy/variance target function and is thus rejecting the update" << std::endl << std::endl;
    vf->from_vector(&init_vars.at(0));
    next_shift = *std::max_element(_h_shifts.begin(), _h_shifts.end());
  }

  if (formic::mpi::rank() == 0)
    formic::of << boost::format("suggested next shift is %.6e") % next_shift << std::endl << std::endl;

  //// load the set of wavefunction variables that gave the lowest energy and choose that shift for the next iteration
  //{
  //  double lowest = shift_improvements.at(0);
  //  int low_i = 0;
  //  for (int i = 1; i < shift_improvements.size(); i++) {
  //    if (shift_improvements.at(i) < lowest) {
  //      low_i = i;
  //      lowest = shift_improvements.at(i);
  //    }
  //  }
  //  if ( lowest > init_energy ) {
  //    rejected = true;
  //    if (formic::mpi::rank() == 0)
  //      formic::of << "linear method is rejecting update because energy increased" << std::endl << std::endl;
  //    vf->from_vector(&init_vars.at(0));
  //    next_shift = *std::max_element(_h_shifts.begin(), _h_shifts.end()) * 10.0;
  //  } else {
  //    vf->from_vector(&shift_wfns.at(low_i).at(0));
  //    next_shift = _h_shifts.at(low_i);
  //  }
  //  if (formic::mpi::rank() == 0)
  //    formic::of << boost::format("suggested next shift is %.6e") % _h_shifts.at(3) << std::endl << std::endl;
  //}

  // reset the shift vector and remember the suggested next shift
  _h_shifts.resize(4);
  _h_shifts.at(3) = next_shift;
  formic::mpi::bcast(_h_shifts);

  return rejected;

}

template bool formic::fqmc::LinMethodUpdater<double>::update(const formic::InputBase & userinp,
                                                             boost::shared_ptr<formic::fqmc::Sampler<double> > & sampler,
                                                             const boost::shared_ptr<formic::CoeffObj<double> > & vf);
template bool formic::fqmc::LinMethodUpdater<std::complex<double> >::update(const formic::InputBase & userinp,
                                                                            boost::shared_ptr<formic::fqmc::Sampler<std::complex<double> > > & sampler,
                                                                            const boost::shared_ptr<formic::CoeffObj<std::complex<double> > > & vf);

/////////////////////////////////////////////////////////////////////////////////////////////////////
///// \brief  Optimize the variables of the supplied coefficient value object via the linear method.
/////
///// \param[in]      userinp  the user's input options
///// \param[in]      co       a smart pointer to the coefficient value object to optimize
/////
/////////////////////////////////////////////////////////////////////////////////////////////////////
//template<class S> void formic::fqmc::linear_method(const formic::InputBase & userinp,
//                                                   const boost::shared_ptr<formic::CoeffObj<S> > & co) {
//
//
//  // print a greeting
//  if (formic::mpi::rank() == 0)
//    formic::of << "Starting linear method optimization." << std::endl << std::endl;
//
//  // print the initial wavefunction variables
//  if ( formic::mpi::rank() == 0 )
//    co->print_variables();
//
//  // get a vector to hold Hamiltonian shifts, as well as the initial shift to use
//  std::vector<double> h_shifts(4, 0.0);
//  h_shifts.at(3) = userinp.get<double>("lm_h_diag_shift");
//
//  // get a vector to hold the energies obtained for different shifts
//  std::vector<double> shift_energies(h_shifts.size(), 0.0);
//
//  // get a vector to hold the wave functions obtained for different shifts
//  std::vector<std::vector<S> > shift_wfns(h_shifts.size());
//  for (int i = 0; i < shift_wfns.size(); i++)
//    shift_wfns.at(i).assign(co->n_der(), formic::zero(S()));
//
//  // initialize the object that will perform the davidson algorithm
//  formic::fqmc::DavidsonLM<S> dsolver(userinp, co);
//
//  // create a vector to hold wavefunction variables
//  std::vector<S> saved_vars(co->n_der(), formic::zero(S()));
//
//  // create a vector to hold the lowest eigenvector in the first derivative subspace
//  std::vector<S> evec(1+co->n_der(), formic::zero(S()));
//
//  // create a list to hold previous updates for use in seeding new krylov subspaces
//  std::list<std::vector<S> > prev_updates;
//  const int max_prev_updates = userinp.get<int>("lm_n_prev_updates");
//
//  // iteratively update the wavefunction
//  for (int iter = 0; iter < userinp.get<int>("max_macro_loop"); iter++) {
//
//    // save a clone of the wavefunction
//    boost::shared_ptr<formic::CoeffObj<S> > old_co = co->clone();
//
//    // save a copy of the wavefunction's variables
//    co->to_vector(&saved_vars.at(0));
//
//    // create an object to hold the starting configuration
//    formic::Configuration start_config(userinp);
//
//    // create an object to hold the starting state of the random number generator
//    formic::LaggedFibonacci<> start_rng;
//
//    // burn in the markov chain
//    if ( !userinp.get<bool>("exact_sampling") ) {
//
//      // prepare a walker
//      formic::VMC_Walker<S> walker(userinp);
//      walker.gf() = co;
//
//      // burn in
//      formic::sample(userinp, userinp.get<int>("lm_burn_len"), walker);
//
//      // save the configuration and the state of the random number generator
//      start_config = walker.config();
//      start_rng = walker.lfg();
//
//    }
//
//    // compute energy and standard deviation of the initial wavefunction for this iteration
//    if (formic::mpi::rank() == 0)
//      formic::of << boost::format("computing energy of initial wavefunction...") << std::endl << std::endl;
//    double init_sdev = 0.0;
//    const double init_energy = formic::fqmc::energy(userinp, 0, userinp.get<int>("lm_samp_len"), start_config, start_rng, co, co, init_sdev);
//
//    // set the shifts that will be tried
//    h_shifts.at(0) = h_shifts.at(3) * 10.0;
//    h_shifts.at(1) = h_shifts.at(3);
//    h_shifts.at(2) = h_shifts.at(3) / 10.0;
//
//    // set the initial guess equal to the initial wavefunction
//    std::fill(evec.begin(), evec.end(), formic::zero(S()));
//    evec.at(0) = formic::unity(S());
//
//    if ( ! userinp.get<bool>("lm_restart_krylov") ) {
//
//      // reset the davidson solver (this also empties the krylov subspace)
//      dsolver.reset(userinp, co, start_config, start_rng, h_shifts.at(0));
//
//      // add the initial guess to the krylov subspace
//      dsolver.add_krylov_vector(userinp, evec);
//
//      // add previous updates to the krylov subspace
//      for (typename std::list<std::vector<S> >::iterator it = prev_updates.begin(); it != prev_updates.end(); it++)
//        dsolver.add_krylov_vector(userinp, *it);
//
//    }
//
//    // compute a wavefunction update and its energy for each of the shifts
//    for (int si = 0; si < h_shifts.size(); si++) {
//
//      // print which shift is being tried
//      if (formic::mpi::rank() == 0)
//        formic::of << boost::format("trying a shift of %.4e") % h_shifts.at(si) << std::endl << std::endl;
//
//      // restore the wavefunction
//      co->from_vector(&saved_vars.at(0));
//
//      // if requested, restart the krylov subspace at each shift
//      if ( userinp.get<bool>("lm_restart_krylov") || dsolver.min_S_eval() < userinp.get<double>("lm_min_S_eval") ) {
//        dsolver.reset(userinp, co, start_config, start_rng, h_shifts.at(0));
//        std::fill(evec.begin(), evec.end(), formic::zero(S()));
//        evec.at(0) = formic::unity(S());
//        dsolver.add_krylov_vector(userinp, evec);
//      }
//
//      // update the krylov subspace for the new shift
//      dsolver.update_ham_shift(userinp, h_shifts.at(si));
//
//      // compute the lowest eigenvector using the davidson algorithm
//      dsolver.solve(userinp, shift_energies.at(si), evec);
//
//      // "normalize" the update
//      for (int i = co->n_der(); i >= 0; i--)
//        evec[i] /= evec[0];
//      //std::vector<S> S_times_evec(co->n_der(), formic::zero(S()));
//      //s_op.operate(co->n_der(), &evec.at(0), &S_times_evec.at(0), formic::mpi::rank(), 0);
//      //for (int i = co->n_der() - 1; i > 0; i--)
//      //  evec[i] /= S_times_evec[0];
//      formic::mpi::bcast(&evec.at(0), evec.size());
//
//      // perform a line search along the update direction
//      formic::fqmc::line_search(userinp, 0, userinp.get<int>("lm_samp_len"), start_config, start_rng, &evec.at(1), userinp.get<double>("lm_max_sdev"), old_co, co);
//
////      // perform a line search along the update direction
////      {
////
////        std::vector<double> step_sizes(1, 0.0);
////        std::vector<double> step_energies(1, init_energy);
////        double best_step_size = 0.0;
////        double best_step_energy = init_energy;
////        int best_step = 0;
////        int step_i = 0;
////        double sdev = 0.0;
////
////        for (int snum = 0; snum < userinp.get<int>("lm_line_search_steps"); snum++) {
////
////          double new_step_size = -1.0;
////
////          // choose what to do next
////          if ( best_step > 0 && best_step < step_sizes.size()-1 && sdev > userinp.get<double>("lm_max_sdev") ) {
////            if (formic::mpi::rank() == 0)
////              formic::of << boost::format("stopping line search due to unexpected too-large standard deviation of %.2e") % sdev << std::endl << std::endl;
////            break;
////          } else if ( sdev > userinp.get<double>("lm_max_sdev") ) {
////            new_step_size = 0.5 * step_sizes.at(step_i);
////            if (formic::mpi::rank() == 0)
////              formic::of << boost::format("erasing most recent step that had a too-large standard deviation of %.2e") % sdev << std::endl << std::endl;
////            step_sizes.erase(step_sizes.begin()+step_i);
////            step_energies.erase(step_energies.begin()+step_i);
////          } else if ( step_sizes.size() < 2 ) {
////            new_step_size = 1.0;
////            if (formic::mpi::rank() == 0)
////              formic::of << boost::format("trying first step...") << std::endl << std::endl;
////          } else if ( step_sizes.size() < 3 ) {
////            new_step_size = step_sizes.at(step_i) * ( best_step_energy == init_energy ? 0.5 : 1.3 );
////            if (formic::mpi::rank() == 0)
////              formic::of << boost::format("trying second step...") << std::endl << std::endl;
////          } else if ( best_step == 0 ) {
////            new_step_size = step_sizes.at(1) * 0.5;
////            if (formic::mpi::rank() == 0)
////              formic::of << boost::format("so far all steps are too long, so we will try a shorter step") << std::endl << std::endl;
////          } else if ( best_step == step_sizes.size()-1 ) {
////            new_step_size = best_step_size * 1.3;
////            if (formic::mpi::rank() == 0)
////              formic::of << boost::format("best step so far was longest step, so we will try an even longer step") << std::endl << std::endl;
////          } else {
////            new_step_size = formic::porabola_min_max(&step_sizes.at(best_step-1), &step_energies.at(best_step-1));
////            if (formic::mpi::rank() == 0)
////              formic::of << boost::format("refining the local minima found during the line search") << std::endl << std::endl;
////          //} else if ( step_energies.at(best_step-1) < step_energies.at(best_step+1) ) {
////          //  new_step_size = 0.5 * ( step_sizes.at(best_step-1) + step_sizes.at(best_step) );
////          //  if (formic::mpi::rank() == 0)
////          //    formic::of << boost::format("refining the local minima found during the line search") << std::endl << std::endl;
////          //} else {
////          //  new_step_size = 0.5 * ( step_sizes.at(best_step+1) + step_sizes.at(best_step) );
////          //  if (formic::mpi::rank() == 0)
////          //    formic::of << boost::format("refining the local minima found during the line search") << std::endl << std::endl;
////          }
////
////          // prepare the next step
////          for (step_i = 0; step_i < step_sizes.size(); step_i++) {
////            if ( new_step_size < step_sizes.at(step_i) )
////              break;
////          }
////          step_sizes.insert(step_sizes.begin()+step_i, new_step_size);
////          step_energies.insert(step_energies.begin()+step_i, 1.0e100);
////          if ( step_i <= best_step )
////            best_step++;
////          //formic::of << boost::format("here 1") << std::endl << std::endl;
////
////          // compute new energy and standard deviation using correlated sampling
////          formic::xscal(co->n_der(), formic::unity(S()) * step_sizes.at(step_i), &evec.at(1), 1);
////          co->update_variables( &evec.at(1) );
////          if (formic::mpi::rank() == 0)
////            formic::of << boost::format("computing energy of updated wavefunction for a %.4e shift and %12.6f step size...") % h_shifts.at(si) % step_sizes.at(step_i) << std::endl << std::endl;
////          step_energies.at(step_i) = formic::fqmc::energy(userinp, 0, userinp.get<int>("lm_samp_len"), start_config, start_rng, old_co, co, sdev);
////          co->from_vector(&saved_vars.at(0));
////          formic::xscal(co->n_der(), formic::unity(S()) / step_sizes.at(step_i), &evec.at(1), 1);
////
////          // if this was the best step size so far, save it
////          if ( sdev < userinp.get<double>("lm_max_sdev") && step_energies.at(step_i) < best_step_energy ) {
////            best_step = step_i;
////            best_step_size = step_sizes.at(step_i);
////            best_step_energy = step_energies.at(step_i);
////          }
////
////          // if the step sizes have become too close together, stop searching
////          const double step_thresh = 1.0e-3;
////          if ( best_step < step_sizes.size()-1 ) {
////            if ( std::abs( step_sizes.at(best_step) - step_sizes.at(best_step+1) ) < step_thresh )
////              break;
////          }
////          if ( best_step > 0 ) {
////            if ( std::abs( step_sizes.at(best_step) - step_sizes.at(best_step-1) ) < step_thresh )
////              break;
////          }
////
////        }
////
////        // if a zero step was best, set the next smallest step as the best
////        if ( best_step == 0 ) {
////            best_step = 1;
////            best_step_size = step_sizes.at(1);
////            best_step_energy = step_energies.at(1);
////        }
////
////        // print result of line search
////        if (formic::mpi::rank() == 0)
////          formic::of << boost::format("best step size was %12.6f with an energy of %20.12f") % best_step_size % best_step_energy << std::endl << std::endl;
////
////        // prepare the update using the best step size
////        formic::xscal(co->n_der(), formic::unity(S()) * best_step_size, &evec.at(1), 1);
////        formic::mpi::bcast(&evec.at(0), evec.size());
////
////      }
////
////      // update the wavefunction's variables
////      co->update_variables( &evec.at(0) + 1 );
//
//      // save this shift's variables
//      co->to_vector(&shift_wfns.at(si).at(0));
//
//      // compute new energy and standard deviation using correlated sampling
//      if (formic::mpi::rank() == 0)
//        formic::of << boost::format("computing energy of updated wavefunction for the %.4e shift...") % h_shifts.at(si) << std::endl << std::endl;
//      double shift_sdev = 1.0e100;
//      shift_energies.at(si) = formic::fqmc::energy(userinp, 0, userinp.get<int>("lm_samp_len"), start_config, start_rng, old_co, co, shift_sdev);
//
//      // if the standard deviation was too large, penalize the energy
//      if ( formic::mpi::rank() == 0 && shift_sdev > userinp.get<double>("lm_max_sdev") ) {
//        shift_energies.at(si) = init_energy + 0.1 * std::abs(init_energy);
//        for (int i = 0; i < si; i++)
//          shift_energies.at(si) += 0.1 * std::abs(shift_energies.at(i));
//        formic::of << boost::format("standard deviation of %.2e is above the user-defined lm_max_sdev limit of %.2e, so the energy is being penalized")
//                      % shift_sdev
//                      % userinp.get<double>("lm_max_sdev")
//                   << std::endl;
//        formic::of << boost::format("penalized energy = %20.12f") % shift_energies.at(si) << std::endl << std::endl;
//      }
//      formic::mpi::bcast(shift_energies.at(si));
//
//      // choose the best shift, using a quadratic fit to interpolate if the interior shift was the best
//      if ( si == 2 && formic::mpi::rank() == 0 ) {
//        const double energy_thresh = 1.0e-9;
//        if        ( shift_energies.at(0) < shift_energies.at(1) && shift_energies.at(0) < shift_energies.at(2) ) {
//          h_shifts.at(3) = h_shifts.at(0) * 3.0;
//        } else if ( shift_energies.at(2) < shift_energies.at(0) && shift_energies.at(2) < shift_energies.at(1) ) {
//          h_shifts.at(3) = h_shifts.at(2) / 3.0;
//        } else if ( std::abs( shift_energies.at(0) - shift_energies.at(1) ) < energy_thresh && std::abs( shift_energies.at(1) - shift_energies.at(2) ) < energy_thresh ) {
//          h_shifts.at(3) = h_shifts.at(1);
//        } else {
//          //const int n = 3;
//          //std::vector<double> lin_reg_mat(n*n, 0.0);
//          //std::vector<double> lin_reg_vec(n, 0.0);
//          //std::vector<double> lin_reg_sol(n, 0.0);
//          //for (int i = 0; i < n; i++) {
//          //  lin_reg_mat.at(i*n+0) = std::log(h_shifts.at(i)) * std::log(h_shifts.at(i));
//          //  lin_reg_mat.at(i*n+1) = std::log(h_shifts.at(i));
//          //  lin_reg_mat.at(i*n+2) = 1.0;
//          //  lin_reg_vec.at(i) = shift_energies.at(i);
//          //}
//          //double d = 0.0;
//          //std::vector<double> work(n*n, 0.0);
//          //std::vector<int> iwork(2*n, 0);
//          //formic::matrix_inverse_lu(n, d, &lin_reg_mat.at(0), &work.at(0), &iwork.at(0));
//          //formic::xgemm('T', 'N', n, 1, n, 1.0, &lin_reg_mat.at(0), n, &lin_reg_vec.at(0), n, 0.0, &lin_reg_sol.at(0), n);
//          //h_shifts.at(3) = std::exp( -0.5 * lin_reg_sol.at(1) / lin_reg_sol.at(0) );
//          std::vector<double> log_shifts(3, 0.0);
//          for (int i = 0; i < log_shifts.size(); i++)
//            log_shifts.at(i) = std::log(h_shifts.at(i));
//          h_shifts.at(3) = std::exp( formic::porabola_min_max(&log_shifts.at(0), &shift_energies.at(0)) );
//          if ( h_shifts.at(3) > 1.1 * h_shifts.at(0) )
//            throw formic::Exception("unexpectedly large shift of %.4e found in the linear regression during the linear method") % h_shifts.at(0);
//          if ( h_shifts.at(3) < 0.9 * h_shifts.at(2) )
//            throw formic::Exception("unexpectedly small shift of %.4e found in the linear regression during the linear method") % h_shifts.at(0);
//        }
//        h_shifts.at(3) = std::max(h_shifts.at(3), userinp.get<double>("lm_min_h_diag_shift"));
//      }
//
//      // broadcast the estimated best shift to all processes
//      if ( si == 2 )
//        formic::mpi::bcast(h_shifts.at(3));
//
//    }
//
//    // save the update for use in seeding future krylov subspaces
//    if (formic::mpi::rank() == 0)
//      prev_updates.push_back(evec);
//    else
//      prev_updates.push_back(std::vector<S>(1));
//    prev_updates.rbegin()->at(0) = formic::zero(S());
//    if ( prev_updates.size() > max_prev_updates )
//      prev_updates.pop_front();
//
//    // load the set of wavefunction variables that gave the lowest energy and choose that shift for the next iteration
//    {
//      double lowest = shift_energies.at(0);
//      int low_i = 0;
//      for (int i = 1; i < shift_energies.size(); i++) {
//        if (shift_energies.at(i) < lowest) {
//          low_i = i;
//          lowest = shift_energies.at(i);
//        }
//      }
//      if ( lowest > init_energy ) {
//        if (formic::mpi::rank() == 0)
//          formic::of << "rejecting update because energy increased" << std::endl << std::endl;
//        co->from_vector(&saved_vars.at(0));
//        h_shifts.at(3) = h_shifts.at(0);
//      } else {
//        co->from_vector(&shift_wfns.at(low_i).at(0));
//        h_shifts.at(3) = h_shifts.at(low_i);
//      }
//    }
//
//    // print the new wavefunction variables
//    if ( formic::mpi::rank() == 0 )
//      co->print_variables();
//
//  }
//
////  // compute final energy
////  if (formic::mpi::rank() == 0)
////    formic::of << boost::format("computing energy of final wavefunction...") << std::endl << std::endl;
////  formic::fqmc::energy(userinp, co);
//
//}
//
//template void formic::fqmc::linear_method(const formic::InputBase &, const boost::shared_ptr<formic::CoeffObj<double> > &);
//template void formic::fqmc::linear_method(const formic::InputBase &, const boost::shared_ptr<formic::CoeffObj<std::complex<double> > > &);
//
/////////////////////////////////////////////////////////////////////////////////////////////////////
///// \brief  Optimize the variables of the supplied coefficient value object via the linear method.
/////
///// \param[in]      userinp  the user's input options
///// \param[in]      co       a smart pointer to the coefficient value object to optimize
/////
/////////////////////////////////////////////////////////////////////////////////////////////////////
//template<class S> void formic::fqmc::linear_method_old(const formic::InputBase & userinp,
//                                                       const boost::shared_ptr<formic::CoeffObj<S> > & co) {
//
//
//  //{
//  //  formic::of << "printing pair compounding:" << std::endl;
//  //  std::vector<int> ioff;
//  //  formic::get_pair_ioff(userinp.get<int>("nao"), ioff);
//  //  for (int j = 1; j < userinp.get<int>("nao"); j++)
//  //  for (int i = 0; i < j; i++)
//  //    formic::of << boost::format("%3i  %3i    %3i") % i % j % ( i + ioff.at(j) ) << std::endl;
//  //    //formic::of << boost::format("%3i  %3i    %3i") % i % j % ( i + userinp.get<std::vector<int> >("orb_pair_ioff").at(j) ) << std::endl;
//  //  formic::of << std::endl;
//  //}
//
//  // print a greeting
//  if (formic::mpi::rank() == 0)
//    formic::of << "Starting linear method optimization." << std::endl << std::endl;
//
//  // print the initial wavefunction variables
//  if ( formic::mpi::rank() == 0 )
//    co->print_variables();
//
//  // get a vector to hold Hamiltonian shifts, as well as the initial shift to use
//  std::vector<double> h_shifts(3, 0.0);
//  h_shifts.at(0) = userinp.get<double>("lm_h_diag_shift");
//
//  // get a vector to hold the energies obtained for different shifts
//  std::vector<double> shift_energies(3, 0.0);
//
//  // iteratively update the wavefunction
//  for (int iter = 0; iter < userinp.get<int>("max_macro_loop"); iter++) {
//
//    // save a clone of the wavefunction
//    boost::shared_ptr<formic::CoeffObj<S> > old_co = co->clone();
//
//    // save a copy of the wavefunction's variables
//    std::vector<S> saved_vars(co->n_der(), formic::zero(S()));
//    co->to_vector(&saved_vars.at(0));
//
//    // create an object to hold the starting configuration
//    formic::Configuration start_config(userinp);
//
//    // create an object to hold the starting state of the random number generator
//    formic::LaggedFibonacci<> start_rng;
//
//    // burn in the markov chain
//    if ( !userinp.get<bool>("exact_sampling") ) {
//
//      // prepare a walker
//      formic::VMC_Walker<S> walker(userinp);
//      walker.gf() = co;
//
//      // burn in
//      formic::sample(userinp, userinp.get<int>("lm_burn_len"), walker);
//
//      // save the configuration and the state of the random number generator
//      start_config = walker.config();
//      start_rng = walker.lfg();
//
//    }
//
//    // compute energy of the initial wavefunction for this iteration
//    if (formic::mpi::rank() == 0)
//      formic::of << boost::format("computing energy of initial wavefunction...") << std::endl << std::endl;
//    const double init_energy = formic::fqmc::energy(userinp, 0, userinp.get<int>("lm_samp_len"), start_config, start_rng, co, co);
//
//    // create a vector to hold the lowest eigenvector in the first derivative subspace
//    std::vector<S> evec(1+co->n_der(), formic::zero(S()));
//
//    // if requested, search for the best shift to use
//    if ( userinp.get<bool>("lm_search_shifts") ) {
//
//      // set the shifts that will be tried
//      h_shifts.at(1) = h_shifts.at(0) * 10.0;
//      h_shifts.at(2) = h_shifts.at(0) / 10.0;
//
//      // compute a wavefunction update and its energy for each of the shifts
//      for (int si = 0; si < h_shifts.size(); si++) {
//
//        // print which shift is being tried
//        if (formic::mpi::rank() == 0)
//          formic::of << boost::format("trying a shift of %.4e") % h_shifts.at(si) << std::endl << std::endl;
//
//        // set the initial guess equal to the current wavefunction
//        std::fill(evec.begin(), evec.end(), formic::zero(S()));
//        evec.at(0) = formic::unity(S());
//
//        // create the Hamiltonian operator
//        formic::fqmc::HamOpLM<S> h_op(userinp, co, start_config, start_rng, h_shifts.at(si));
//
//        // create the overlap matrix operator
//        formic::fqmc::OvlpOpLM<S> s_op(userinp, co, start_config, start_rng);
//
//        // create the preconditioning operator
//        formic::fqmc::PrecOpLM<S> p_op(userinp, co, start_config, start_rng, shift_energies.at(si), h_shifts.at(si));
//
//        // compute the lowest eigenvector using the davidson algorithm
//        formic::start_timer("davidson");
//        formic::davidson(1+co->n_der(),
//                         userinp.get<int   >("lm_davidson_iter"),
//                         userinp.get<double>("lm_davidson_thresh"),
//                         shift_energies.at(si),
//                         &evec.at(0),
//                         &h_op,
//                         &s_op,
//                         &p_op,
//                         false);
//        formic::stop_timer("davidson");
//
//        // update the wavefunction's variables
//        for (int i = co->n_der(); i >= 0; i--)
//          evec[i] /= evec[0];
//        //std::vector<S> S_times_evec(co->n_der(), formic::zero(S()));
//        //s_op.operate(co->n_der(), &evec.at(0), &S_times_evec.at(0), formic::mpi::rank(), 0);
//        //for (int i = co->n_der() - 1; i > 0; i--)
//        //  evec[i] /= S_times_evec[0];
//        formic::mpi::bcast(&evec.at(0), evec.size());
//        co->update_variables( &evec.at(0) + 1 );
//
//        // compute new energy using correlated sampling
//        if (formic::mpi::rank() == 0)
//          formic::of << boost::format("computing energy of updated wavefunction for the %.4e shift...") % h_shifts.at(si) << std::endl << std::endl;
//        shift_energies.at(si) = formic::fqmc::energy(userinp, 0, userinp.get<int>("lm_samp_len"), start_config, start_rng, old_co, co);
//
//        // restore the wavefunction
//        co->from_vector(&saved_vars.at(0));
//
//      }
//
//      // choose the best shift, using a quadratic fit to interpolate if the interior shift was the best
//      if (formic::mpi::rank() == 0) {
//        if        ( shift_energies.at(1) < shift_energies.at(0) && shift_energies.at(1) < shift_energies.at(2) ) {
//          h_shifts.at(0) = h_shifts.at(1);
//        } else if ( shift_energies.at(2) < shift_energies.at(0) && shift_energies.at(2) < shift_energies.at(1) ) {
//          h_shifts.at(0) = h_shifts.at(2);
//        } else {
//          const int n = 3;
//          std::vector<double> lin_reg_mat(n*n, 0.0);
//          std::vector<double> lin_reg_vec(n, 0.0);
//          std::vector<double> lin_reg_sol(n, 0.0);
//          for (int i = 0; i < n; i++) {
//            lin_reg_mat.at(i*n+0) = std::log(h_shifts.at(i)) * std::log(h_shifts.at(i));
//            lin_reg_mat.at(i*n+1) = std::log(h_shifts.at(i));
//            lin_reg_mat.at(i*n+2) = 1.0;
//            lin_reg_vec.at(i) = shift_energies.at(i);
//          }
//          double d = 0.0;
//          std::vector<double> work(n*n, 0.0);
//          std::vector<int> iwork(2*n, 0);
//          formic::matrix_inverse_lu(n, d, &lin_reg_mat.at(0), &work.at(0), &iwork.at(0));
//          formic::xgemm('T', 'N', n, 1, n, 1.0, &lin_reg_mat.at(0), n, &lin_reg_vec.at(0), n, 0.0, &lin_reg_sol.at(0), n);
//          h_shifts.at(0) = std::exp( -0.5 * lin_reg_sol.at(1) / lin_reg_sol.at(0) );
//          if ( h_shifts.at(0) > h_shifts.at(1) )
//            throw formic::Exception("unexpectedly large shift of %.4e found in the linear regression during the linear method") % h_shifts.at(0);
//          if ( h_shifts.at(0) < h_shifts.at(2) )
//            throw formic::Exception("unexpectedly small shift of %.4e found in the linear regression during the linear method") % h_shifts.at(0);
//        }
//      }
//      formic::mpi::bcast(h_shifts.at(0));
//      h_shifts.at(0) = std::max(h_shifts.at(0), userinp.get<double>("lm_min_h_diag_shift"));
//
//    }
//
//    // print the shift that we chose
//    if (formic::mpi::rank() == 0)
//      formic::of << boost::format("estimate for best shift is %.4e") % h_shifts.at(0) << std::endl << std::endl;
//
//    // compute and apply the update using the chosen shift
//    {
//
//      // set the initial guess equal to the current wavefunction
//      std::fill(evec.begin(), evec.end(), formic::zero(S()));
//      evec.at(0) = formic::unity(S());
//
//      // create the Hamiltonian operator
//      formic::fqmc::HamOpLM<S> h_op(userinp, co, start_config, start_rng, h_shifts.at(0));
//
//      // create the overlap matrix operator
//      formic::fqmc::OvlpOpLM<S> s_op(userinp, co, start_config, start_rng);
//
//      // create the preconditioning operator
//      formic::fqmc::PrecOpLM<S> p_op(userinp, co, start_config, start_rng, shift_energies.at(0), h_shifts.at(0));
//
//      // compute the lowest eigenvector using the davidson algorithm
//      formic::davidson(1+co->n_der(),
//                       userinp.get<int   >("lm_davidson_iter"),
//                       userinp.get<double>("lm_davidson_thresh"),
//                       shift_energies.at(0),
//                       &evec.at(0),
//                       &h_op,
//                       &s_op,
//                       &p_op,
//                       false);
//
//      // update the wavefunction's variables
//      for (int i = co->n_der(); i >= 0; i--)
//        evec[i] /= evec[0];
//      //std::vector<S> S_times_evec(co->n_der(), formic::zero(S()));
//      //s_op.operate(co->n_der(), &evec.at(0), &S_times_evec.at(0), formic::mpi::rank(), 0);
//      //for (int i = co->n_der() - 1; i > 0; i--)
//      //  evec[i] /= S_times_evec[0];
//      formic::mpi::bcast(&evec.at(0), evec.size());
//      co->update_variables( &evec.at(0) + 1 );
//
//      // compute new energy using correlated sampling
//      if (formic::mpi::rank() == 0)
//        formic::of << "computing the new wavefunction's energy..." << std::endl << std::endl;
//      shift_energies.at(0) = formic::fqmc::energy(userinp, 0, userinp.get<int>("lm_samp_len"), start_config, start_rng, old_co, co);
//
//      // if the energy went up, undo the update
//      if ( shift_energies.at(0) > init_energy ) {
//        if (formic::mpi::rank() == 0)
//          formic::of << "rejecting update because energy increased" << std::endl << std::endl;
//        co->from_vector(&saved_vars.at(0));
//      }
//
//    }
//
//    // print the new wavefunction variables
//    if ( formic::mpi::rank() == 0 )
//      co->print_variables();
//
//  }
//
////  // compute final energy
////  if (formic::mpi::rank() == 0)
////    formic::of << boost::format("computing energy of final wavefunction...") << std::endl << std::endl;
////  formic::fqmc::energy(userinp, co);
//
//}
//
//template void formic::fqmc::linear_method_old(const formic::InputBase &, const boost::shared_ptr<formic::CoeffObj<double> > &);
//template void formic::fqmc::linear_method_old(const formic::InputBase &, const boost::shared_ptr<formic::CoeffObj<std::complex<double> > > &);
